Phase behavior of hard-core lattice gases: A Fundamental Measure approach 
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We use an extension of fundamental measure theory to lattice hard-core fluids to study the phase 
diagram of two different systems. First, two-dimensional parallel hard squares with edge-length 
a — 2 in a simple square lattice. This system is equivalent to the lattice gas with first and second 
neighbor exclusion in the same lattice, and has the peculiarity that its close packing is degenerated 
(the system orders in sliding columns). A comparison with other theories is discussed. Second, 
a three-dimensional binary mixture of parallel hard cubes with ox = 6 and as = 2. Previous 
simulations of this model only focused on fluid phases. Thanks to the simplicity introduced by 
the discrete nature of the lattice we have been able to map out the complete phase diagram (both 
uniform and nonuniform phases) through a free minimization of the free energy functional, so the 
structure of the ordered phases is obtained as a result. A zoo of entropy-driven phase transitions 
is found: one-, two- and three-dimensional positional ordering, as well as fluid-ordered phase and 
solid-solid demixings. 



I. INTRODUCTION 

Hard-core systems are the paradigm of entropy-driven 
phase transitions. The first example of an entropy-driven 
(orientational) ordering transition is given in the famous 
Onsager's papesi on the isotropic-nematic transition in 
a three-dimensional system of thin hard rods. But prob- 
ably, the best known and discussed example of entropy- 
driven (three-dimensional positional) ordering transition 
is the freezing of hard spheres. This was first devised by 
Kirkwood et al^ using an approximate theory, but the 
definite evidence about the existence of such a purely en- 
tropic transition was the numerical simulations of Alder 
and Wainwright, 3 and Wood and Jacobsoni^ As very few 
models can be solved exactly, definite conclusions on the 
existence of phase transitions often come from numeri- 
cal simulations. But in many cases these are very de- 
manding and powerful computers are needed in order to 
reach a reliable system size. This fact, together with the 
inexistence of appropriate theoretical approaches, could 
explain that until the end of the eighties there were no 
more instances of entropy-driven ordering transitions. At 
that time, a series of numerical simulations 5 ^^ 7 " showed 
that hard-core interaction can also induce one- and two- 
dimensional positional ordering (smectic and columnar 
phases, respectively, in liquid crystal terminology). This 
was a very striking fact, because it was generally believed 
that the mechanism underlying these phase transitions 
was the decrease of internal energy rather than the gain 
of entropy. 

Apart from ordering transitions, it is well known that 
binary nonadditive mixtures can demix by a pure en- 
tropic effect. An extreme case of nonadditivity was stud- 
ied by Widom and Rowlinson 8 in a model with two dif- 
ferent species interacting ideally between members of the 
same species (ctaa = obb = 0) and with a hard-core in- 
teraction between unlike particles (ctab = o~). They rig- 
orously showed that the system demixes into two fluid 



phases with different compositions. This can be eas- 
ily understood if we notice that the available volume is 
more effectively filled by pure phases than by the mix- 
ture. Another interesting example of this kind is found 
in colloid-polymer mixtures. Experimentally, it is well 
known that the addition of non-adsorbing polymers to 
a colloidal suspension induces an effective attraction be- 
tween the colloidal particles that can induce the floccula- 
tion of the colloid. A simple explanation for this effect is 
that the clustering of colloids (large particles) leaves more 
free volume to the polymers (small particles) , what trans- 
lates into a gain of entropy. This mechanism is known as 
depletion. Many models2*ifl*i 1 « have been successfully in- 
troduced in order to illustrate how this effect can induce 
a fluid-fluid phase separation in mixtures. 

Special mention merits the case of the additive binary 
mixture of hard spheres. The absence of a spinodal in- 
stability in the Percus-Yevick solution for this systemic 
led to believe that entropic demixing was not possible 
for additive mixtures. But almost thirty years later, 
Biben and HasenA 3 . predicted such a spinodal by using a 
more accurate integral equation theory. Since this result, 
a bunch of theoretical, 14 i 15 . 16 i 17 i 18 i 19 simulation 20 ! 21 ! 22 ! 23 
and experimentalS 4 *^ 5 ^ 6 *^! results appeared supporting 
the existence of demixing in additive binary mixtures of 
hard spheres when the diameter ratio is at least 5:1. Al- 
most at the same time, it was pointed out that instead of 
a fluid- fluid demixing at least one of the separated phases 
might be ordered, 16 ! 17 ! 18 ! 19 ! 21 ! 22 ! 23 ! 25 ! 26 ! 27 The actual see- 
nario for this system is a mctastablc fluid-fluid demixing 
(not confirmed by direct simulation 22 ) preempted by a 
fluid-solid coexistence or (if the mixture is sufficiently 
asymmetric) a solid-solid one. Qualitatively, this is the 
same situation one finds in a binary additive mixture of 
parallel hard cubesi 28 i 29 i 30 

From a theoretical point of view one of the first exactly 
solvable hard-core models showing a fluid-solid transi- 
tion was a lattice model proposed by Temperley^ Many 
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other lattice hard-core models were studied in the six- 
ties by adapting the approximate theories developed for 
Ising-like models to hard body systems. 32 ^ They suc- 
ceeded in the prediction of an order-disorder transition 
and mainly focused on studying the dependence of the 
nature of the transition upon the range of the hard-core 
and the topology of the underlying lattice. 

For the continuum model of hard spheres, one of the 
most successful theories to study the freezing has been 
density functional theory (DFT). Many accurate func- 
tionals have been devised for the monocomponcnt fluid, 34 
but when they are applied to the binary mixture some 
problems arise: (i) many of the theories employed are 
not directly formulated for mixtures and the extension is 
far from being straightforward, (ii) it is very difficult to 
study the solid phases because it is not trivial to deter- 
mine which is the most stable structure for the mixture, 
and this information is an input in most approaches^ 
These difficulties have been circumvented by mapping the 
binary hard-core mixture into a monocomponent fluid 
(large particles alone) with a hard-core and an effec- 
tive short-range attractive potential. It is then possible 
to use perturbation theory in order to study the phase 
diagram^ On the other hand, the solid phase is usu- 
ally assumed to be an fee crystal of the large particles 
with the small particles uniformly distributed. Although 
this approach has been extensively ugedj^iiSiSSiSLiSi it is 
only valid for low molar fractions of the small particles. 
Besides, even in this case, the assumption that the den- 
sity of small particles is uniform in the ordered phase is 
rather unrealistic because the ordering of large particles 
induces structure in the distribution of the small ones. 
To the best of our knowledge, this problem has not been 
addressed satisfactorily yet. 

A direct study (without resorting to an effective one- 
component fluid) has been carried out for a binary mix- 
ture of parallel hard cubes with Rosenfeld's fundamental 
measure theory 28 i 29 i 30 (FMT). This theory has the ad- 
vantage of being naturally formulated for mixtures. A 
complete analysis of fluid-fluid demixing has been per- 
formed for arbitrary size ratios, but again, the lack of 
intuition about the distribution of small particles in the 
crystal makes impossible to study freezing in this sys- 
tem. To solve this problem one should perform a free 
minimization of the free energy functional and obtain the 
structure of the ordered phases as an output. But due to 
the continuum nature of the system this would require a 
huge amount of numerical work. 

The situation is more favorable for the lattice counter- 
part of this model. Indeed, simulations of a binary mix- 
ture of parallel hard cubes (6 : 2) on a simple cubic lat- 
tice were performed by Dijkstra and Frenkel^ but their 
focus was whether entropic demixing could be observed 
in additive binary mixtures and the structure of the in- 
homogeneous phases was not considered. The results of 
these simulations (a stable fluid-fluid demixing) are in 
contradiction with the predictions of the continuum sys- 
tem (it exhibits a fluid-fluid spinodal only for size ratios 



above 10 : 1, and it is always preempted by freezing of the 
large component). With the aim of explaining this mis- 
match, we extended the fundamental measure functional 
for parallel hard cubes to the lattice version^ With this 
theory we have shown in a previous work^i that the latter 
is the correct picture. (A more detailed account of this 
work will be given here.) Furthermore, due to the dis- 
crete nature of the system it is possible to give a complete 
description of the ordered phases (see below). Thus, we 
show that lattice models, treated in a suitable manner, 
can serve as a starting point to study the structure of 
ordered phases in (continuum) mixtures. 

There is a second benefit of this extension of FMT 
to lattices that we want to emphasize. These simula- 
tions, together with an exactly solvable model proposed 
by Widorc 9 and Frenkel and Louis 10 in different contexts, 
show that lattice models can give accurate descriptions 
of demixing phenomena. But in spite of their historical 
role in the development of Statistical Mechanics and their 
simplicity, with a few exceptions density func- 

tional theories have only focused on continuum models. 
We believe that the formulation of classical density func- 
tional approaches for lattice models will help to better 
understand both, the phase behavior of complex fluids 
and the formal structure of the approximate functionals. 

The paper is organized as follows. A review of the lat- 
tice version of FMT is presented in Sec.[n] In Sec. lIIII we 
use this theory to obtain the complete phase diagram of 
two systems. First, a two-dimensional system of parallel 
hard squares with edge-length a = 2 on a square lattice 
(this is equivalent to the two-dimensional lattice gas with 
first and second neighbor exclusion). This system has 
been widely studied in the literature (see Refs. |33 and 
l33l and references therein) and there exists a big con- 
troversy about its phase behavior so far unsettled. A 
detailed analysis is performed by applying the new the- 
ory, and a comparison with results from other theories 
is discussed. The lattice fundamental measure theory 
(LFMT) appears to be at the same level of accuracy of 
the other well accepted theories. Secondly, we have ad- 
dressed the problem of the binary additive mixture stud- 
ied by simulations, i.e. a binary mixture of parallel hard 
cubes (<tl = 2, <jg = 6) on a simple cubic lattice. Due to 
the simplification introduced by the lattice, the complete 
phase diagram has been mapped out. It shows a very 
rich collection of entropic phase transitions. As a matter 
of fact, we have found one-, two- and three-dimensional 
ordering transitions, as well as fluid-ordered phase and 
solid-solid demixings. A free minimization of the free 
energy functional has been performed, so the structure 
of the ordered phases has also been obtained. Finally, 
conclusions are discussed in Sec. II VI 



II. THEORY 

The construction of the lattice fundamental measure 
functional is based on the ideas of the exact zero- 
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dimensional reductionist and in the exact form of the 
one-dimensional functional. A full account of the details 
of the procedure can be found in Ref. ^3- In that work, 
the general form of the functional for a system of parallel 
hard cubes in a simple cubic lattice (for any dimension, 
particle size or number of components) is presented. Ba- 
sically, the idea behind it is to construct a family of func- 
tionals for arbitrary dimension in such a way that they 
consistently satisfy the dimensional reduction property of 
the exact functionals down to zero-dimensional cavities 
(i.e. cavities which can host no more than one particle). 
Moreover, the prescription chosen is inspired in the ex- 
act functional for the one-dimensional system, which is 
recovered with the scheme proposed. 

Let us consider a d-dimensional additive mixture of 
parallel hypercubes with edge- lengths +e lattice 

spacings, where a is the species index and e = 0, 1 does 
not depend on a, i.e all the species have, simultaneously, 
even or odd sizes (the mixed — nonadditive — case is more 
involved^ so we just ignored it because it will not be 
used anywhere in this work) . In Ref. the excess free 
energy functional for this system in this approximation 
was found to be [cf. Eq. (3.2) of that reference] 



sGZ d ke{0,l} d 



d-k 



o(« (k) ( S )), (1) 



where k — &/i ®o(v) = V + (1 — 77) ln(l — rj) is 

the excess free energy for a zero-dimensional cavity with 
mean occupancy O77 < 1, the reciprocal temperature 
in Boltzman's units and n^(s) are weighted densities 
defined by the convolutions 



(2) 



p a (s) being the one-particle distribution function for 
species a and 



1=1 



1 if — a a — k — e < s < a D 
otherwise. 



(3) 
(4) 



Notice that as weights are indexed by k € {0, there 
are 2 d different weighted densities. 

The direct correlation function between species a and 
7 can be obtained from this functional as 
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where n k = J2 a a> ak G a — l) d k Pa are the weighted densi- 



ties Q in the uniform limit and tp^y (s) is the convolution 



(6) 



re: 



Because of the structure of the direct correlation function 
it is convenient to work with its discrete Fourier trans- 
form, which takes the form 

e« 7 (q) = - E y^^H-qX'fq), (7) 

ke{o,i} d k 



where 
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sin qi/ 2 



(8) 



The general expression of the functional adopts 
very simple forms when particularized to specific sys- 
tems. In order to make clear the structure of the func- 
tional, we will introduce a diagrammatic notation which 
helps visualizing its dimensional reduction properties in 
a simple way. For the sake of simplicity let us consider 
the lattice gas with first and second neighbor exclusion 
in a two-dimensional square lattice. This is a system of 
parallel hard squares with a = 2 lattice spacings. In dia- 
grammatic notation the excess free energy functional l(T)l 
can be written 

(3F^[p] = E[*o(n)- <i> o(o^)-$o(^)+$o(o)], 

sez 2 

(9) 

where the diagrams represent the weighted densities @ 
as 



,(i,i) 



0^3 =n (1 < 0) (s) 



o = 



(s) = p(s 1 ,s 2 ) + p{si + l,s 2 ) 

+p{si, s 2 + 1) + p(si + 1, s 2 + 1), 
p(si, s 2 ) + p(si + 1, s 2 ), (10) 
(s) = p(si,s 2 ) +p(si,s 2 + 1), 
(s) = p(si,s 2 ). 



What becomes apparent with this diagrammatic no- 
tation is that the excess functional Q can be regarded 
as a linear combination of contributions due to a partic- 
ular set of zero-dimensional cavities I|1U|) . Furthermore, 
we can manipulate the diagrams in order to prove the 
dimensional reduction properties that the functional @ 
satisfies. To illustrate this, we will consider the dimen- 
sional reduction to a one-dimensional system, the hard 
rod lattice gas, whose exact excess functional is known 
to have the form (Q. 40 To perform this reduction, we 
will apply an infinite external potential in every site of 
the lattice except in an infinite linear chain defined by 
C = {(si,0) : si £ Z}. This implies that the centers of 
mass of the particles can only occupy the sites in £, the 
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system becoming equivalent to a hard rod lattice gas with 
particles of size a — 2. In terms of p(s) this means that 
p(s) = p(si)5 S2j o, where Sij is the Kronecker symbol and 
p(s) is the one-particle distribution function for the one- 
dimensional system. Within this constraint, the excess 
free-energy functional of the effective system can be ob- 
tained by summing over S2 £ Z in (O- Each contribution 
gives, respectively, 



»2t 



E $ o(n) = 2$ o (<~), 

s 2 ez 

®o{^<>) = $ (o^), 

= 2$ o(o), 

- 2 GZ 

^ $o (o) = $0(0), 



s 2 ei 



where the diagrams in the r.h.s. must now be interpreted 
as 0-0 = p(s) + p(s + 1) and o = p(s). Therefore, 
the excess free-energy functional for the one-dimensional 
system so obtained is 

P^l d) \p\= £[*o(<~>)-*o (<>)], (11) 



which coincides with the exact result [see Eqs. (2.26) and 
(3.1) of Ref.EJ. 

Another example in three-dimensions is given in the 
Appendix. 



III. APPLICATIONS 

In spite of the simple structure of the lattice funda- 
mental measure functionals, the applications to specific 
systems have proven able to describe very complex phase 
diagrams^! In this section we will study in detail two 
particular systems: first, the lattice gas with exclusion 
to first and second neighbors on a square lattice; second, 
a binary mixture of parallel hard cubes with a = 2 and 
6, in a simple cubic lattice. 



A. Parallel hard squares 

This model is defined by the interaction pair potential 



00 if \si - < 1 for both i = 1, 2, 
otherwise. 



(12) 



It has been previously studied employing other approxi- 
mate theories, such as finite-size scaling methodsf 4 ^^ se- 
ries expansions22i3L£i and clusters methods All au- 
thors agree in that the close-packed state is a columnar 
phase (ordered along one dimension but fluid along the 



other). This notwithstanding, the nature of the transi- 
tion remains doubtful, the results being highly depen- 
dent on the theory usedi^ While some authors conclude 
that the system exhibits a third order transition very 
near close packing, 47 - 49 j£2i£i others obtain a second or- 
der transition at a lower densityj^MLS and even some 
of them have speculated about the lack of such a phase 
change The results obtained with the present the- 
ory are in accordance with those of the second order 
phase transition. Unfortunately, we have no concluding 
arguments to umpire this dispute. 

The advantage of our approach over other theories is 
that it provides a simple prescription to build a density 
functional in closed form. Then, all the powerful tools 
of density functional theory may be applied. The excess 
functional for this system within the present theory is 
that of Eq. ©. 

In order to pin down the phase diagram of the sys- 
tem we will proceed systematically: first, studying the 
uniform phases, and second, considering spatial inhomo- 
geneities. For a uniform density profile p(s) = p, the 
weighted densities H1U|) become = Ap, 00 = £ = 2p 
and o = p, where < p < 1/4. The excess free energy 
density (in k^T units) can be calculated particularizing 
©, which yields 



$ ox (p) = $o(4p) - 2$ (2p) + $ (p), 



(13) 



Adding up the ideal gas contribution, <j>id = p(ln p — 1), 
and taking into account the definition of $0, we obtain 
the following free energy density for the fluid 



$fluid(p) =plnp+ (1 - p)ln(l - p) 
+ (1 - Ap) ln(l - Ap) 
- 2(1 - 2p) ln(l - 2p). 



(14) 



From this, all the thermodynamic properties of the fluid 
phase can be derived. For instance, the fugacity is given 
by 



Zfluid 



p(l - 2pY 



(l-p)(l-4p)*' 
and the pressure takes the simple form 

(l-2p) 2 



/?Pfluid = hi 



(l-p)(l-4p) 



(15) 



(16) 



The structure of the equilibrium phase can be ana- 
lyzed by means of the direct correlation function, ob- 
tained from JSJ particularizing for d — 2 and a single 
component with a = 2. In a symmetry broken continu- 
ous phase transition, for some q ^ 



1 - pc(q) = 0, 



(17) 



this condition being equivalent to the divergence of the 
structure factor. Since we are interested in the spatial 
instabilities of the uniform phase, we have to look for the 
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FIG. 1: The unit cell for a columnar phase with periodicity 
equal to two lattice spacings is shown in (a), and for a solid 
phase with the same periodicity in (b). 



lowest value of p which makes the condition 1|17J) solvable 
for some q. Taking into account that the symmetry of the 
system enables us to take q = (q, 0), Eq. (fTTjl becomes 



cos 2 (g/2) 



4p j P_ 

1— 2p ' 1-p 



-ip 



l-2p 



(18) 



Since the denominator is positive in the whole range of 
p, and so is the numerator for small values of the density, 
the solution corresponds to the vanishing of the latter. 
This occurs at the density 



Pcrit 



3- V5 



0.1910, 



(19) 



and, of course, q = n, implying the periodicity of the in- 
homogeneous phase to be d — 2ir/q = 2 lattice spacings. 

What remains to be determined is the symmetry of 
the nonuniform phase at the transition point. Based on 
the previous results for this system and on those recently 
obtained for parallel hard cubes in the continuum^ we 
guess that this phase must be either a columnar or a 
solid (ordering along the two coordinate axes). In order 
to determine which phase is the stable one, we have per- 
formed a global minimization of the functional (0 within 
the constraints imposed by the symmetry and periodic- 
ity of both the columnar and the solid phases. For a 
generic columnar phase with periodicity equal to two lat- 
tice spacings, the one-particle distribution function takes 
the form 



Pcoi(s) 



pi if si is even, 
p2 otherwise, 



(20) 




FIG. 2: Sublattice densities for the columnar phase. It is 
also shown the metastable fluid beyond the transition point. 



takes the form 



$coi (pi ; p) = $id {pi , 2p - pi ) + $ (4p) - $o (2p) 



\\Mtpi 

$oG°i) - 



) + $o(4p-2p!) 
*o(2p-pi)], 



(22) 



where we have substituted the density profile (|20|l in JJjJ , 
used the relation 2p = p\ + p2, and introduced the ideal 
term $id(pi,P2) = ^Y,iPiO- n Pi - l )- 

We can now minimize the total free energy density at 
constant p. Note that in this case we have to minimize 
with respect to a single variable. The Euler-Lagrange 
equation is 



p 1 (l-2 Pl ) 2 (l-2p + p 1 ) 
(2p-pi)(l-4p + 2pi) 2 (l-pi) 



= 1. 



(23) 



One solution corresponds to the uniform phase (p^ q = p) . 
It can be easily checked that this is indeed the minimum 
of the free energy for p < p cr ;t, as expected. After re- 
moving this solution from l|23|l . we obtain a quadratic 
polynomial whose roots become physical for p > p cr it- 
Above the critical density the uniform phase is no longer 
a minimum; instead, we have a columnar phase given by 



while for a solid phase with the same periodicity we have 



pi if s\ and S2 are even, 
Psol(s) = { p2 if si or s 2 is odd, 
P3 otherwise. 



(21) 



A sketch of the unit cell for each case is shown in Fig. ^ 
Note that the uniform phase is included in both l|20|l and 

For the columnar phase the total free energy density 



1 /(l-2p)(p-p crit )(3-2p-2p crit ) 



eq 

* =P+ 2 



(we have chosen p^ q > p!j q — 



(24) 



2p— pj q ). This phase has a 
lower free energy than the fluid phase for p > p cr ;t, but 
we still have to calculate free energy for the solid branch 
in order to know which one is the stable phase above the 
transition point. 

For the solid phase, substituting l|21|l in (0 and adding 
the ideal contribution, the total free energy density turns 
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out to be 

^soi (pi , Ps ; p) = $id (pi , 2p - p+ , p 3 ) + $0 (4p) 
-$ (2p + p_)-$ (2p-p-) 

+ - [* (pi) + 2$ (2p - p+) + $o(p 3 ) 



(25) 



where we have used 4p = p\ + 2p2 + P3 to eliminate the 
dependence on p2, and have defined p± — (pi ± P3)/2. 
As in the previous case, the equilibrium density profile 
is the global minimum of i|25|) at constant p, but now 
we have two independent variables, pi and ps. Thus, the 
Euler-Lagrange equations are now the system of algebraic 
equations 



Pl (l-2p-p_) 2 (l-2p + p + ) 

(i-p 1 )(2p- P+ )(i-2p + p_y 

p 3 (l-2p + p_) 2 (l-2p + p + ) 
(l-p 3 )(2p-p+)(l-2p-p_) 2 



(26) 



= 1. 



The fluid phase, given by p° q = pg q = p, is the solution 
for p < p cr it- The solution for p > p CTlt must be obtained 
numerically and is plotted in Fig. 01 In Fig. 0] we can see 
that the solid branch bifurcates with a free energy lower 
than that of the fluid phase, but larger than that of the 
columnar phase. The transition is then fluid-columnar. 

It is feasible to study analytically the behavior of each 
branch at the transition point. This would give a definite 
conclusion about the nature of the phase change. It is 
straightforward to check the continuity of d$/dp at p cr i t 
(for both the solid and the columnar branches), but a 
discontinuity is found in the second derivative at p cr i t , so 
the transition is second order. Furthermore, the stable 
phase beyond the transition point is the one with lowest 



second derivative for p 



$ fluid(piit) 



^colVf crit/ 
^sollPcrit) 



■ Pcrif From the values 

2(15 + 7^5) -61.3, 
2(5 + 2\/5) w 18.9, 
4(5 + V5) w 28.9, 



(27) 



we conclude that indeed the system undergoes a second 
order transition from a fluid phase to a columnar one 
at Pcrit- Besides, as it can be inferred from the density 
dependence of the free energy density for every branch 
(Fig. 01, the columnar phase remains the most stable 
phase up to close packing. 

We can now plot the equation of state (Fig. |SJ, with 
the fluid branch given by Eq. (|16fl and the columnar one 
by 



(l-2p) 2 (l-2^)(l-4p + 2p c 1 q ) 
(l-4p)2(l-^)(l-2p + p e ^) 



The fugacity of the columnar phase is given by 

_ ^(l-2p) 2 (l -2p^) 2 
Zco1 (1 - 4p) 4 (l 

At the critical point, we have 



eq\ 

Pi) 



PPcrit = hi 2, 



11 + 5V5 



(28) 



(29) 



(30) 



As mentioned at the beginning of this section, the re- 
sults from this lattice fundamental measure theory are 
compatible with the ones obtained by Bellemans and 
Nigamia (p crit » 0.202, /3p crit « 0.788 and z cl . it « 17.29) 
through the cluster method of Rushbrooke and Scoins 54 
(plotted with a dashed line in Fig. EJ. From Fig. we 
can see that the agreement at low and high densities is 
very accurate, and deviations occur only near the critical 
point. This can be understood if we realize that both the- 
ories neglect correlations beyond a certain distance be- 
tween the particles, so a very accurate description of the 
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FIG. 5: Equation of state of the hard square fluid from lattice 
fundamental measure theory (solid line) and from the cluster 
method of Rushbrooke and Scoins (dashed line). 



critical properties should not be expected. This notwith- 
standing, as remarked by Runnels^ due to the degener- 
acy of the close-packed configuration, this system is diffi- 
cult to study with finite-size or series expansions method, 
and a closed-form approximation could be superior at de- 
scribing the correct phase behavior. 



B. Multicomponent System of Parallel Hard Cubes 

Let us now consider a multicomponent hard cube lat- 
tice gas in a simple cubic lattice. If we denote a±, . . . , a p 
the edge-lengths of the different species, then the inter- 
action potential between species a and 7, will be given 
by 



6 a ^y (S, S ) 



00 if max Is* 

i=l,2,3 

otherwise. 



(31) 

The lattice fundamental measure approximation for the 
free energy functional of this system has already been 
reported in Ref. |4fJ [Eqs. (3.2) and (3.3) of that refer- 
ence] , together with the phase diagram for the particular 
case of a binary mixture with ctl = 6 and as = 2, but 
no details about the calculation were given. In this sub- 
section we will study the phase behavior of the general 
uniform mixture and obtain the complete bulk phase di- 
agram, including both uniform and ordered phases, for 
the particular case just mentioned. 

In the uniform regime, the one-particle distribution 
function no longer depends on the spatial variables: 
p a (s) = p a , a = 1, . . . ,p. In this case, the free energy 



$(pi, ...,p p ) = ^2p a {\np a - 1) + $o(n 3 ) 

a=l 

-$o(n 2 ) + $o(rci)-$o(no) (32) 

[with the densities n% defined below Eq. JSJ]. 

As it is well known, the stability of the mixture is de- 
termined by the matrix 



'ay 



<9 2 $ 

dp a dp 1 



(33) 



In order for the system to be stable in a mixed state this 
matrix must be positive definite. As this requirement is 
fulfilled in the low density limit, the spinodal curve can be 
determined through the equation det M = 0. When the 
excess free energy density of the system depends on the 
densities only via the finite set of moments £j = J2 a a aPa 
(I = 0, ...,771), the spinodal can be expressed in the 
equivalent, but more suitable formS det Q = 0, where 



Q, 



4\ 



fe=0 



d 2 <P cx 



(34) 



This is just our case, because the excess free energy den- 
sity [the $0 contributions in depends on the densi- 
ties through the set {n^}, and this variables can easily 
be expressed in terms of the set of moments {£o> • • • , £3}- 
Thus the equation for the spinodal of our system reads 



(1 + 2&) 2 - (6 + 3$0(1 + 2&) 

-3(&-6)(l + £4)+6(5& 



£0=0. (35) 



For a binary mixture with the small component having 
as = 2, it can be shown that the smallest size ratio, 
r = at,/ as, necessary to have a spinodal instability is r = 
13. This value is in strong disagreement with previous 
simulation results^ which reported a demixing phase 
transition for r — 3. An explanation of this mismatch 
has already been provided in Ref. l4ll and will become 
clear later on. 

Some spinodals for different size ratios are shown 
in Fig. It should be noticed that the continuum 
counterparts^ is recovered in the limit as — * 00, while 
keeping r constant. For that system, it was shown that 
the minimum value of r to find demixing is r = 5 + v24 fts 
9.98. Therefore, we can conclude that the lattice en- 
hances the stability of the mixture. What this analogy 
with the continuum system suggests is that we should ex- 
pect fluid-fluid demixing to be preempted by the freezing 
of one of the coexisting phases also in the lattice model. 
It must be remarked that, unlike in the continuum case, 
in the lattice system the stability condition involves not 
only the size ratio, but also the edge-length of one of the 
components, thus making the analysis of the stability 
more complex. 
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FIG. 6: Spinodal curves for a uniform binary mixture with 
the smallest component of size erg = 2 and different values of 
the size ratio, r — <tl/<ts- 

The discrete nature of this system provides a very suit- 
able framework to study ordering transitions. In the con- 
tinuum, only a partial analysis have been done^S because 
the minimization of the functional becomes a numeri- 
cally very demanding task. Usually this minimization is 
performed by restricting the density profiles into a para- 
metric class, and then minimizing with respect to one 
or a few parameters. The problem for the mixture is 
that it is very difficult to guess the appropriate class. In 
contrast, the situation in the lattice is easier to handle 
because it is feasible to perform a free minimization with 
the only constraints imposed by the symmetry of the or- 
dered phase and its periodicity. As it was shown in the 
analysis of the two-dimensional system, the periodicity of 
the ordered phase can be estimated from the divergence 
of the structure factor, in the case of a mixture the latter 
being a matrix. The analog to the condition (II Tf) for the 
mixture is 

det(P" 1 - C(q)) = 0, (36) 

where P = (6 a ^p a j, a diagonal matrix, and C(q) = 
(c« 7 (q)) is the matrix of Fourier transforms of the di- 
rect correlation functions between all species. 

In the remaining of this subsection we will restrict our- 
selves to the particular case of the binary mixture with 
erg = 2 and ctl = 6, the only case for which simula- 
tions are available^ The main result of these simula- 
tions is that the mixture undergoes an entropy-driven 
fluid-fluid demixing, thus being the only known example 
of an athermal additive model showing this feature. 

The strategy we have adopted in order to obtain the 
complete phase diagram for the mixture has been the 
following: (i) First, we have calculated the phase di- 
agram for the pure component systems, both for the 
small and large particles; (ii) then, we have obtained the 
curves marking spatial instabilities for the whole mix- 
ture through the divergence of the structure factor ma- 
trix, and have calculated the periodicity of the ordered 
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FIG. 7: Equation of state (pressure vs. packing fraction) for a 
cr = 2 parallel hard cube lattice gas on a simple cubic lattice. 
The different symmetries are denoted by F, Sm and C, mean- 
ing fluid, smectic and columnar, respectively. The periodicity 
has been indicated by a subindex. The inset shows a very 
narrow first order transition from a smectic to a columnar 
phase. 

phases arising at the bifurcation points, and (iii) finally, 
we have completed the phase diagram by calculating all 
the possible phase transitions compatible with the results 
obtained in the two previous steps, choosing those ther- 
modynamically more stable. 

For the pure component systems, we will proceed as in 
the two-dimensional case. For a = 2, the fundamen- 
tal measure excess free energy functional Q is given 
in diagrammatic notation in l|A.I|) . For a general one- 
component system with edge-length cr, the functional 
form is obtained particularizing but the structure 
is the same as that of (H3J). 

As it was mentioned in the previous subsection, in a 
symmetry broken continuous phase transition condition 
(|T7Jl is satisfied. For cr = 2 it yields 

?7crit = cr 3 Pcrit = 0.568, (37) 

with q = tt indicating a periodicity d = 2n/q = 2, while 
for cr = 6 we find the value 

Vent = <J 3 Pcrit = 0.402, (38) 

and period d — 7. Notice that at close packing this sys- 
tem has period d = 6, so we must consider both. With 
respect to the symmetry of the phases, we have consid- 
ered the smectic, columnar and solid ones (ordering along 
one, two and three coordinate axes, respectively). 

Now it is possible to perform a free minimization of 
the functional with the above restrictions. In this case 
we have to proceed numerically, because the complex- 
ity of the problem does not permit an analytical treat- 
ment. This notwithstanding, the structural form of the 
functional simplifies the numerical work: the weighted 
densities are just convolutions which can be computed 
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FIG. 8: Equation of state (pressure vs. packing fraction) for a 
(7 = 6 parallel hard cube lattice gas on a simple cubic lattice. 
The different symmetries are denoted by F, Sm, C and S, 
meaning fluid, smectic, columnar and solid, respectively. The 
periodicity has been indicated by a subindex. 
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FIG. 9: Free energy density of the fluid (dotted line), smectic 
(dashed line), columnar (solid line) and solid (dotted-dashed 
line) phases, for the system of hard cubes with a = 2 in a 
simple cubic lattice. 



by using fast Fourier transform. To give an idea about 
the degree of complexity of the problem we will say that 
the simplest phase to minimize is the period-2 smectic, 
which involves minimization on two variables, and the 
most complex one is the period-7 solid, which involves 
minimization on twenty variables. 

The phase diagrams of both systems, a — 2 and a = 6, 
are shown in Figs. and [SJ respectively. Also, the free 
energy density near the critical point appears in Figs. El 
and ^3 respectively. From these figures, we can see 
that there exists a strong competition between the dif- 
ferent phases. This reflects in the very narrow first order 
transitions observed in the phase diagrams, such as the 
Sni2-C2 coexistence (J7sm 2 — 0.673 and r/c 2 — 0.677) 
in the a = 2 system, and the S7-C7 (i]s 7 = 0.617 and 
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Tl 

FIG. 10: Free energy density of the fluid (dotted line), smec- 
tic (double dotted-dashed line), columnar (dotted-dashed 
line) and solid (double dashed-dotted line) phases with peri- 
odicity d = 7 and smectic (long-dashed line), columnar (solid 
line) and solid (dashed line) phases with periodicity d = 6, 
for the same system of the previous figure but with a = 6. 



i]c 7 — 0.631) for cr = 6. Since our treatment is approxi- 
mate, these phase transitions could actually be spurious: 
given the small differences between the free energy den- 
sities of the phases involved, other scenarios might be 
possible. In contrast, there also exists very well defined 
transitions which offer higher confidence, such as the C7- 
C6 first order transition (ric 7 — 0.656 and r]c 6 — 0.827) 
in the a = 6 system. 

Let us now consider the binary mixture with <tl = 6 
and as = 2. From the discussion about the uniform 
multicomponent system we can conclude that, for a size 
ratio r = 3 and a small particle edge- length as — 2, there 
is no fluid- fluid demixing (not even metastable). Then, 
we have to look for spatial instabilities. To this purpose 
we must study condition (|36|l . The direct correlation 
function is now a 2 x 2 matrix whose elements are given 
by (jsj). If we characterize the thermodynamics of our 
system by the total packing fraction rj = r/^ + i]s, where 
^(S) = a L(S)P L (S) i s ^ ne packing fraction of the large 
(small) cubes, and by the composition x — rj-^/rj, then, 
for every value of x, we have to look for the smallest 
value of rj which makes the condition (|36|l solvable. The 
solution is plotted in Fig. ^] (Fig. ^| shows the same 
phase diagram with rj replaced by the pressure) . It should 
be remarked that for < x < 0.728 the period of the 
ordered phases at the fluid spinodal is d = 2, while for 
1 > x > 0.728 we have found d = 7. 

With this guidance we can start looking for the true 
phase diagram. This is a very demanding numerical task, 
but feasible in a reasonable amount of time. For each 
coexistence curve involving, say, phases Pi and P2, we 
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FIG. 11: Phase diagram of the binary mixture of parallel 
hard cubes (size ratio 6 : 2), r\ = t?l + r/s being the total 
packing fraction of the large (L) and small (S) cubes, and 
x = T%lr\. The phases are labeled F (fluid), Sm a (smectic), 
C a (columnar), and S a (solid), where a = 2,6,7 stands for 
the periodicity of the ordered phases. The dashed lines join 
coexisting states. The dotted line corresponds to the spin- 
odal of the uniform fluid. For 0.81 < x it marks a stable 
continuous F-S7 phase transition. The circles are coexisting 
states taken from the simulation results in Ref. Finally, 
the stars correspond to the states whose density profiles are 
represented in Figs. 113^161 

have to solve the equilibrium equations 

Zl{Vp 1 ,x Fi ) = 2 L (??p 2 ,2;p 2 ), (39) 
Zs07Pi,£Pi) = ^sOyps; 2 ^)- 

Every iteration of the procedure requires the minimiza- 
tion of the functional, at constant rj and x, for both 
phases. In the simplest case this corresponds to a min- 
imization problem with four variables, but in the most 
complicated case we have to deal with a forty variable 
minimization. Another problem we find is that the subtle 
differences between the free energy of different phases, al- 
ready encountered in the monocomponent systems, make 
it very hard to discern which one is the most stable phase. 
So, in many cases it is not clear which coexistence is 
thermodynamically more stable. In these doubtful cases, 
we have resorted to the Gibbs free energy per particle, 
g(p, X) = X[1l + (1 — X)fi§, where /xl,s is the chemical 
potential of each species and X = p^/ ' p = xj (r+(l — r)x) 
is the molar fraction. When g(p, X) is plotted at con- 
stant pressure as a function of X, the coexisting phases 
in the mixture can be found through a double tangent 
construction. 

The complete phase diagram appears in Figs. ^JJ and 
n"21 in two different representations. In Fig. it can 

be observed (see the insets) that there exist very small 
parts on the phase diagram with a plethora of very nar- 
row coexistence regions. As in the monocomponent case, 



many of them could just be spurious. One of the most 
remarkable features is that there exists a wide phase 
separation between a small-particle-rich fluid phase and 
a large-particle-rich columnar phase (which becomes a 
solid phase for higher pressures). As explained in the 
Introduction, this is the usual scenario for this kind of 
mixtures. The revision of the simulations (also shown 
in Fig. 1 1 1|) resulting from this phase behavior has al- 
ready been discussed in detail in Ref. \^[ The main con- 
sequence of this comparison is that some of the state 
points obtained in the simulations must have been mis- 
interpreted as a fluid, while they should exhibit columnar 
ordering. Another interesting result is that, at approxi- 
mately (3p « 0.24, there appear a solid-solid phase sepa- 
ration between a small-particle-rich and a large-particle- 
rich phases. Finally, it is worth mentioning the existence 
of an extremely narrow chimney of Sg-Cg coexistence. 
This suggests that the columnar phase is very sensitive 
to small perturbations. 

As concerns the density profile in the ordered phases, 
we have chosen a few representative states of the system 
(marked with stars in Fig. Illfl in order to illustrate the 
way large and small particles distribute in each phase. In 
Fig. El we show the density profiles for a period-7 solid 
phase for different sections perpendicular to the direction 
S3. This size is one lattice spacing bigger than the size 
of the large particles. We can conclude from the figure 
that the large particles are arranged in such a way that 
the small ones can be accommodated between any two 
of them. So, we could say that the large particles in 
the unit cell are surrounded by the small ones. Note 
that the density of the latter does not change very much 
within the unit cell. In Fig. ^] the density profile for 
a period-6 solid is shown. It is very interesting to see 
that its structure is completely different to the previous 
one. Now, as the unit cell is of the same size of the large 
particles, the small ones can only be placed at vacant 
unit cells. Thus, the unit cell is completely filled with 
the small particles, which, as can be appreciated in the 
figure, form a crystal. Noticeably, the value of the small 
particle density is slightly higher at the contact with the 
large particles, which can be interpreted as an adsorption 
phenomenon. Figures [TBI and [TBI exhibit the profile of the 
period-6 column and smectic phases, respectively. From 
the discussion of the period-6 solid profiles we have just 
made, the interpretation of these two new density profiles 
should be straightforward. 



IV. CONCLUSIONS 

In this work we have applied the extension of FMT 
to lattice models^ to study two systems already treated 
in the literature: the two-dimensional lattice gas with 
first and second neighbor exclusion on a square lattice, 
and the binary mixture of parallel hard cubes with edge- 
lengths (Tl = 6 and as = 2 on a simple cubic lattice. 
For both systems we have shown a systematic way of us- 
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ing the theory in order to perform a complete calculation 
of the bulk phase diagrams. In particular, for the first 
system a very detailed analysis have been carried out. 
All the relevant thermodynamic functions can be ana- 
lytically obtained within this approach, for both uniform 
and columnar phases. The results compare well with oth- 
ers previously obtained by different authors using widely 
accepted theoretical approaches. It is remarkable that, 
in spite of the simplicity of the treatment, the results ob- 
tained are rather accurate, specially in the low and high 
density limits. 

On the other hand, a very complete study has been 
carried out for the three-dimensional system. The rich 
phase diagram obtained is very striking, considering the 
simplicity of both the system and the theory. There 
appear many different entropy-driven phase transitions: 
fluid-ordered phase demixing, one-, two- and three- 
dimensional ordering and solid-solid phase separation. 
The results obtained have also allowed to reinterpret the 
simulations results^ in a way consistent with the general 
picture emerged during the last decade about demixing 
of additive hard-core binary mixtures: fluid-fluid demix- 
ing is always preempted by the ordering of one of the 
phases. 

We have performed a free minimization of the func- 
tional and have thus obtained the structure of all the 
ordered phases. The results we get show that the density 
profile of the small particles is far from being uniform. 
In general, when dealing with three-dimensional models, 
a free minimization is not feasible, so the density profile 
must be properly parametrized. Since there is very lit- 
tle intuition about the density profiles of mixtures, the 
small components are always assumed to be uniformly 
distributed over the volume. According to our findings, 
this is definitely wrong. Perhaps our results can help to 



gain insight into what a proper parametrization of the 
density profiles looks like. 

We would also like to emphasize that the LFMT is a 
mean- field- like theory, but not a trivial one. This can be 
understood if we realize that the direct correlation func- 
tions have finite range, which means that at some point 
the correlations between particles are neglected. This 
fact is reflected in the lower accuracy of the description 
around critical points. However, away from these regions 
the results are far more accurate than those obtained 
from a standard mean field theory, and what is even more 
important, at the expese of no much more work. 

The study of lattice gases within the framework of den- 
sity functional theory can be very fruitful, as we hope to 
have been able to transmit in this work. On the one hand, 
there is not loss in phenomenology or complexity of be- 
haviors; on the other hand, the approach is considerably 
simpler numerically, something that allows one to tackle 
problems which have so far not been solved in continuum 
models. Besides, from a purely formal point of view, lat- 
tice density functionals reveal some structures which may 
be hidden in similar developments for continuum models. 
So, its careful analysis may reveal important properties 
of the latter in the near future. For these reasons we 
believe that it is very important to extend this kind of 
work to more general lattice models. We plan to report 
on that shortly. 
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(a) Large particles 





(b) Small particles 



FIG. 13: Density profiles for a solid phase with periodicity 
d = 7, corresponding to r\ — 0.56, x = 0.97 and /3p = 0.02. 
Different planar sections at S3 = 1, 2, 3, 4 (from top to bottom 
and from left to right) are plotted. 
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APPENDIX 



Let us consider a system of parallel hard cubes in a simple cubic lattice with edge- length a = 2. This system is 
equivalent to the lattice gas on a simple cubic lattice with first, second and third neighbor exclusion. The excess 
functional, in diagrammatic notation, can be written as 

PT C M= E [ $0 (c^) - $ o(n)- $ o(^)- $ »(^) + +^(H + *o(^)-$o(o)] , (A.l) 

sez 3 

where the diagrams represent the weighted densities 



n (W) (s)= V p( Sl ,s 2 + i,s a +.j), o-o=nf°' 1 '°) 



i,j,fc=0,l 

(s) = E /°( S lj S 2 + M 3 +j) 

j,i=o,i 

hj=0.1 

J? =n( 1 '°' 1 )(s) = E p( s i+hS2,s 3 +j), 
1,4=0,1 

In Ref. ^3 it was pointed out that the family of approximate functionals constructed after the prescription of the 
lattice fundamental measure theory consistently satisfied the dimensional reduction property of the exact functionals. 
Hence, the functional JHJ must be recovered from the functional (| A. 1|> . To show it we can apply an infinite external 
potential in every site on the tridimensional lattice, except in the sites laying on the plane V = {(0, s 2 , S3) : s 2 , S3 € Z}. 



(s) 




f«), 




i=0,l 




(s) 




S3), 




i=0,l 




(s) 




S3), 




i=0,l 




(s) 


= p(si,s 2 , S3). 
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Then the different contributions to the effective excess functional become 

E* (^H*°(n). e*°(s)=*°(s)> 

siGZ siGZ 

E $ o(o^°) = 24> (o^), ^ $ (^) = 2$ (o), 

siez sigz 

£ $ (jj) = 2 *o(^) , ]T$o(o) = $o(o), 

si£Z siGZ 

and therefore (|A.1(1 reduces to ©. 

This diagrammatic notation can also be extended to multicomponent systems, although for more than two compo- 
nents it becomes too cumbersome. I 



